Effect of quantization of vibrations on the structural properties of crystals 
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We study the structural effects produced by tfie quantization of vibrational degrees of freedom in 
periodic crystals at zero temperature. To this end we introduce a methodology based on mapping a 
suitable subspace of the vibrational manifold and solving the Schrodinger equation in it. A number of 
increasingly accurate approximations ranging from the quasi-harmonic approximation (QHA) to the 
vibrational self-consistent field (VSCF) method and the exact solution are described. A thorough 
analysis of the approximations is presented for model monoatomic and hydrogen-bonded chains, 
and results are presented for a linear HF chain where the potential energy surface is obtained 
via first-principles electronic structure calculations. We focus on quantum nuclear effects on the 
lattice constant, and show that the VSCF is an excellent approximation, meaning that correlation 
between modes is not extremely important. The QHA is excellent for covalently-bonded, mildly 
anharmonic systems, but it fails for hydrogen-bonded ones. In the latter, the zero-point energy 
exhibits a non-analytic behavior at the lattice constant where the H-atoms center, which leads to a 
spurious secondary minimum in the quantum-corrected energy curve. An inexpensive anharmonic 
appoximation of non-interacting modes appears to produce rather good results for hydrogen-bonded 
chains, for small system sizes. However, it converges to the incorrect QHA results for increasing 
size. Isotope effects are studied for the first-principles HF chain. We show how the lattice constant 
and the HF distance increase with decreasing mass, and how the QHA proves to be insufficient to 
reproduce this behavior. 
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I. INTRODUCTION 

Structural properties of solids like interatomic dis- 
tances, bond angles and equilibrium lattice parameters 
are customarily calculated by assuming that atomic nu- 
clei (or ionic cores) behave as classical particles. Within 
this framework, the electronic problem is solved for a 
fixed configuration of clamped nuclei. The resulting 
ground state electronic density allows for the computa- 
tion of the forces acting on the nuclei through Hellman- 
Feynman's theorem. These forces arc then used to deter- 
mine the equilibrium configuration. While this is justified 
for a large variety of systems of interest, the assumption 
of classical nuclei becomes questionable whenever light 
atoms such as hydrogen are involved. Moreover, this ap- 
proximation can not address isotope effects; the ground 
state electronic energy does not depend on the nuclear 
masses but only on the atomic numbers. A widely-used 
approach to introduce the quantum nature of the nuclei 
is the quasi-harmonic approximation (QHA)jii2ii^ In this 
approximation, the nuclear ground state energy is sup- 
plemented with the zero-point-energy (ZPE) correspond- 
ing to harmonic nuclear vibrations. This is based on 
a second-order Taylor expansion of the potential energy 
surface (PES) in terms of the atomic coordinates around 
the equilibrium configuration of the nuclei. The resulting 
dynamical matrix, given by the second derivatives of the 
potential with respect to the nuclear coordinates, can be 
diagonalized thus leading to a set of orthogonal eigenvec- 
tors or normal modes Q , together with the corresponding 
set of frequencies ui* 

To simulate an infinite crystal, it is customary to apply 
periodic boundary conditions on the unit cell.^ Therefore, 



a crystal described by a basis of N atoms in the unit cell is 
characterized by 3A^ vibrational bands ui (k) , where k is 
a wave vector in the phonon Brillouin zone (BZ). Within 
this framework, the QHA energy at zero temperature, 
j^QHA^ is defined as 



3N „ 

^ / gik)u;i{k,V)dk (1) 



where V is the volume of the system, ^(k) is the den- 
sity of states and the integral extends to the Brillouin 
zone. The first term in the RHS is the ground state en- 
ergy for the equilibrium configuration at volume V. The 
second term includes the quantum nature of the nuclei 
through the harmonic ZPE, which also depends on V in 
a way that is characteristic of the type of bonding. For 
covalent bonding it increases upon compression, while 
for hydrogen-bonding there is a competition between in- 
creasing and decreasing frequencies. The combination of 
ground-state and zero-point energies leads to a modified 
]t;QHA Y Q^j-ve - in contrast to E'^^ vs. V, with a 
minimum located at a volume corrected by quantum nu- 
clear effects. The QHA has the additional advantage of 
lending itself naturally to the simultaneous incorporation 
of quantum and thermal effects.^ 

In the harmonic approximation the ground state nu- 
clear wave function is a product of single- mode Gaussians 
on the vibrational normal mode variables Q , and the ex- 
pectation value of the nuclear coordinates coincides ex- 
actly with the classical nuclear configuration. Apart from 
variations mediated by volume changes, internal struc- 
tural parameters (distances and angles) are unaffected in 
the QHA. It is then expected that the QHA breaks down 
when treating highly anharmonic systems. 
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Additional terms in the Taylor expansion lead to an- 
harmonicity, which can come essentially in two forms. 
The simplest one is intra-mode anharmonicity, when the 
vibrational modes remain non-interacting but the poten- 
tial felt by one or more modes cannot be approximated 
by a quadratic expression. This type can be included in 
a rather simple way. Normal modes can be used as coor- 
dinates to map the PES beyond the harmonic level, and 
the energy of each mode obtained by solving a set of in- 
dependent one-dimensional Schrodinger equations. The 
second type of anharmonicity is due to mode-coupling. 
When the excursions along a vibrational normal coor- 
dinate are large, the approximation of small oscillations 
implicit in the second-order expansion breaks down. Suc- 
cessive terms in the expansion involve products of modes 
of the form C™C;" well as higher-order terms involving 
more than two modes. Which groups of modes are more 
strongly coupled depends very much on the specific sys- 
tem. In general there are no rules to find out a priori 
which couplings need to be considered. A possible strat- 
egy is to displace the system along one mode, say Ci, by 
an amount e, and then optimize the atomic coordinates 
under the constraint Ci = £• By projecting the displace- 
ments obtained for the optimized configuration onto the 
original normal modes, one can identify and select the 
modes with the largest projection.^ The main goal of this 
paper is to evaluate the quality of the QHA and, where 
required, to provide improved schemes at an affordable 
computational cost. To this end, we will include anhar- 
monicity by solving the vibrational Schrodinger equation 
in a sequence of increasingly accurate approximations, 
and compare to the QHA results. 

An additional issue that arises when computing the 
ZPE is that the BZ integrals in ([T]) must be replaced 
with appropriate averages over a finite set of representa- 
tive points in the BZ. This can be done in several ways. 
The most accurate one is to compute the force constants 
in real space and then use them to obtain the dynam- 
ical matrix at an arbitrarily dense mesh of k-points in 
the BZ. This requires calculations in large supercells or, 
alternatively, linear response calculations On the other 
extreme, the crudest approximation is to include only the 
BZ-center modes of the unit cell, while improvements can 
be achieved by considering BZ-center modes of larger su- 
percells. The advantage of this latter is that it is easier 
to extend to anharmonic situations. A second goal of 
this paper is to analyze how large a supercell should be 
in order to reproduce the equilibrium structure obtained 
with converged BZ averaging. This issue has rarely been 
discussed in the literature. 

In order to answer these questions we have studied 
two qualitatively different models: a one-dimensional 
monoatomic chain and a diatomic linear chain that mod- 
els a double-well hydrogen-bonded system. In Section 
nil Al we show that, as usually done in the literature, sys- 
tems characterized by covalent bonding can be safely de- 
scribed within the quasi-harmonic approximation using a 
rather coarse BZ sampling. In Section [ill Bl in contrast. 



we show that hydrogen-bonded systems require some 
level of anharmonicity in their treatment. An anomaly 
in the ZPE appears when the protons center in the H- 
bonds, upon compression. This produces an unusual be- 
havior of the QHA, which can lead to a secondary spuri- 
ous minimum in the energy-volume curve. We then study 
a realistic linear F-H- • -F chain where the PES has been 
obtained from first-principles calculations fSection llll C\ . 
We study quantum effects, in particular the isotope ef- 
fect on equilibrium lattice constant and internal geome- 
try (H-F distance) occurring when the nuclear masses are 
modified fSection llll E[) . In this case, an inexpensive an- 
harmonic approximation that neglects mode-coupling ap- 
pears to produce an accurate lattice constant. In Section 
IIVI we present our conclusions and elaborate on possible 
extensions. In the following, we introduce the theoretical 
and computational approaches used in the present work. 



II. METHODS AND APPROXIMATIONS 

In this paper we will only consider one-dimensional 
systems, for the sake of simplicity. Although computa- 
tionally more demanding, extensions to higher dimen- 
sions are straightforward. In the QHA, the energy for a 
one-dimensional system containing A'' particles is 



N-l 



(2) 



where a is the lattice constant - which plays the role 
of the volume, and E''''{a) is the ground state eletronic 
energy in the approximation of classical nuclei (classical 
energy). The vibrational frequencies uji{a) are obtained 
by assuming periodic boundary conditions (PBC), i.e. 
we consider a cyclic chain. The sum runs up to iV — 1 
because the vibrational mode corresponding to the rigid 
translation of the whole chain has zero frequency. 

In the harmonic approximation there is no need to 
consider explicitly the iV-atom chain. The same result 
can be obtained by studying monoatomic or diatomic 
cells subject to PBC. This equivalence is explained thor- 
oughly in most solid state books. ^ Here we just quote 
the main result, i.e. the dispersion relations that express 
the frequency as a function of wavenumber k in the one- 
dimensional phonon Brillouin zone. For a monoatomic 
chain the dispersion relation is 



a) = 2 



V"{a) 



sm 



(3) 



with m the mass of the particles, k — Tm/{N — l)a and 
n — 0,1, ■ ■ ■ , N — 1. The quantity V"{a) is the second 
derivative (curvature) of the interaction potential V eval- 
uated at the equilibrium configuration, and depends on 
the lattice constant a. The frequencies obtained in this 
way correspond exactly with those of an A^-atom cyclic 
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chain. In the hmit of N 
monoatomic chain is 



cxD, the QHA energy for a 



E^"^{a) = E''\a) + ^ j g{k) uj{k, a) dk 
2 Jbz 



(4) 



The case of a diatomic hnear chain is analogous, with the 
difference that there are two dispersion relations, repre- 
senting the acoustic and optical branches. The expres- 
sions are still analytic, but slightly more complicatedi^ 

For a hydrogen-bonded diatomic chain like F-H- • -F 
the harmonic approximation can only be carried out at 
one of the two equivalent global minima. If the barrier 
is so high that tunnelling is unlikely, this is a reasonable 
approximation. However, if tunnelling is important, an- 
harmonic effects start to play a significant role and the 
QHA breaks down. By symmetry, one could also choose 
a reference configuration at the top of the inversion bar- 
rier, i.e. with the H-atoms centered between F-atoms. 
This configuration is a stationary point of the PES, but 
the QHA breaks down because a whole portion of the 
optical phonon branch is unstable around zone-center. 

The calculated normal modes can now be used to 
map the true PES, beyond the harmonic approxima- 
tion. In general, this is a complicated mutidimensional 
fitting problem that requires specific techniques such as 
the product representation^. A much simpler alterna- 
tive is to keep considering the normal coordinates as 
non-interacting, but anharmonic. This corresponds to 
pure intra-mode anharmonicity. The full-dimensional vi- 
brational problem is then reduced to mapping the one- 
dimensional potential for every normal coordinate keep- 
ing the other normal coordinates at their classical value, 
i.e. zero. Then, the 1-D Schrodinger equation is solved 
for each one of these potentials and their ground state 
energies are added up. We call this the anharmonic ap- 
proximation (ANHA). 

At the other end of the theoretical spectrum, the prob- 
lem of interacting phonons can be tackled exactly by solv- 
ing the {N — I)-dimensional Schrodinger equation. This 
task can be accomplished with little effort for three de- 
grees of freedom, but finds a hard wall at six. This means 
that, in the one-dimensional case, one can routinely com- 
pute the exact solution for a 4-atom supercell, but can- 
not go beyond 7-atom chains if a significant fraction of 
the eigenvalue spectrum is required. Stochastic methods 
such as diffusion Monte Carlo can, however, be used to 
compute the low-energy region of the spectrum, as it has 
been done for vibrations in molecular systemsi^ 

Here we have chosen the strategy of solving the vi- 
brational Schrodinger equation using a discrete vari- 
able representation. We used Lagrange grids based 
on a combination of Cartesian and Hermite orthogonal 
polynomialsiSiii for double-well and single-well poten- 
tials, respectively. When solving the Schrodinger equa- 
tion exactly, a maximum of three vibrational degrees of 
freedom was considered. This implies that the largest 
monoatomic linear chain we studied was made of four 
unit cells, while the largest diatomic chain contained two 



unit cells. 

We have also studied a class of approximations to the 
full-dimensional solution inspired by the vibrational self- 
consistent field (VSCF) method,— which is frequently 
used to compute vibrational excitation spectra of large 
molecules In the VSCF method the normal coordinates 
are taken as completely uncorrelated and the total wave 
function is written as a product of single- mode wave func- 
tions: 



*(ci,---,cjv-i)= n '^■'(^■' 



(5) 



in the same spirit of the Hartree approximation to the 
problem of interacting electrons. Contrary to the har- 
monic and anharmonic approximations, each mode now 
feels the presence of the other modes, but in a mean-field 
way. The VSCF equations are 



(6) 



where 



V, 



N-l 

(CO = y • • • y ^(Ci, • • • , Cn-i) n I -^^ (0) P dC: 



(7) 

are the mean-field potentials for each mode. For a I-D 
system the energy in the VSCF method is given by 



VSCF 



N-l 

E 



{N 



v^iC) I MO P dC (8) 



where the second term corrects for double-counting in the 
sum of eigenvalues. Notice that the integral in this term 
should be independent of z, within numerical accuracy. 

The determination of the mean-field potential requires 
the computation of (iV — 2)-dimensional integrals, which 
rapidly becomes an expensive operation. Therefore, it is 
desirable to find either simpler approximations or some 
way of evaluating these integrals at a reduced compu- 
tational cost. The anharmonic approximation (ANHA) 
mentioned above removes this limitation by replacing the 
single- mode densities \(t>j{Cj)\'^ in with delta- functions 
centered at (^j — 0. With this, the ANHA can be carried 
onto much larger systems. 

The quality of the VSCF results depends on the choice 
of coordinates. If the system is only weakly anharmonic, 
then the classical normal modes are a very reasonable 
choice. This is not necessarily the case when anharmonic- 
ity becomes important. The question arises of whether 
it is still possible to find a set of orthogonal coordinates 
that are only weakly correlated, or whether correlation 
is intrinsic to the problem and cannot be significantly 
reduced by a clever choice of coordinates. Within the 
VSCF context, finding the optimally de-correlated co- 
ordinates is equivalent to requiring that these minimize 
EvscF- The optimal coordinates {Ci} can be obtained 
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by solving the VSCF problem for a linear combination 
of normal modes CI = J2j ^ij Cj j ^-nd minimising Eysc f 
with respect to the coefficients Zij of the rotation matrix, 
subject to orthonormality constraints.^'' This task can be 
accomplished by constrained multidimensional optimisa- 
tion algorithmSfi^ Nevertheless, there is no clear evidence 
of a general method for obtaining the optimal set of co- 
ordinates. In this work we analyze the quality of the 
VSCF approximation with respect to the choice of vi- 
brational coordinates for the ab initio hydrogen-bonded 
linear chain. 



III. RESULTS 
A. Model monoatomic chain 

We first considered a periodic monoatomic chain where 
the atoms interact via a Morse potential of the form 



V{x) = oil-e-'' 



(9) 



where parameter D represents the binding energy of a 
dimer, ag is the location of the minimum of the poten- 
tial, which corresponds to the classical lattice constant, 
and 6 is a parameter that determines the curvature of 
the potential through the relation V"{aQ) — 2Db^. The 
Morse potential is intrinsically anharmonic, and allows 
us to study the effect of quantization of the vibrations 
on the lattice constant. For the various approxima- 
tions described in the previous Section we calculated the 
quantum-corrected energy as a function of the lattice pa- 
rameter a, and then determined the minimum of the E 
vs. a curve. 

In the solution of Schrodinger's equation, the extent of 
quantum nuclear effects can be measured by the product 
7 = mD, with m the mass of the particles and D the 
energy scale of the potential. Small values of 7 lead to 
important quantum effects, becoming less relevant as 7 
increases, and eventually converging towards the classical 
results for 7 —>■ 00. Therefore, we report our results 
in terms of this quantity. In this work we have used 
ao = 1.3983 Bohr and b = 2.9864 Bohr-^ 

We first analyzed the effect of Brillouin zone sampling 
on the equilibrium lattice constant. For the QHA this 
can be done for an arbitrarily large number TV of k-points 
in the BZ, which is equivalent to considering a supercell 
containing an equal number TV of unit cells. In Fig. [T] 
we show the convergence pattern of the QHA lattice con- 
stant as a function of 7, for increasing values of TV. 

It is obvious that quantum effects are more important 
for smaller masses and for weaker interaction potentials. 
For H-atoms and D ^ 2 eV (7 — 135 a.u.) we are 
in the region to the left of Fig[TJ where the quantum- 
corrected lattice constant is about 4% larger than the 
classical value (horizontal line). It can be seen that, for 
this type of potential, 10 cells are already sufficient to 




y = mD (a.u.) 

FIG. 1: (Color online) QHA lattice constant as a function of 
7, for increasing number of atoms in the supercell: TV = 2 
(black, thin solid), TV = 3 (red, short-dashed), TV = 4 (green, 
long-dashed), TV = 10 (blue, dot-dashed), and TV = 2000 
(orange, thick solid-light grey). The dashed horizontal line 
represents the classical value. 



reproduce the infinite crystal lattice constant to high ac- 
curacy. However, 4 cells already produce excellent results 
and 2 cells underestimate a in less that 1%. 

In order to illustrate the origin of the lattice expan- 
sion, we show in Fig. [5] the various energy curves for 
a 4-atom supercell, as a function of the lattice parame- 
ter for 7=304.97 a.u. The classical energy curve E'^''{a) 
(black, solid line) exhibits a minimum at the equilib- 
rium lattice constant, where the one-dimensional pres- 
sure, P ~ -dE'^\a)/da, vanishes. The ZPE (red, dot- 
dashed line) depends on the lattice constant through the 
curvature of the potential. Upon compression (reducing 
a), the chain becomes stifFer, all frequencies increase and 
so does the ZPE. Therefore, the QHA energy given by 
the sum of the two terms in (j4|), and represented by the 
blue, dashed curve, exhibits a minimum at an expanded 
lattice constant. This effect is similar to thermal expan- 
sion, but the origin is purely quantum-mechanical. It can 
also be ascribed to a zero-point pressure originated in the 
dependence of the ZPE on the lattice constant. 

To evaluate the quality of the QHA we calculated the 
equilibrium lattice constant in the ANHA and VSCF ap- 
proximations, and exactly for a 4-atom supercell. The 
differences between these results and the QHA lattice 
parameter, Aa, are reported in Fig. [3l It can be seen 
that these differences are one order of magnitude smaller 
than those arising from a poor sampling of the BZ. Thes 
observations justify the approach commonly used to de- 
termine the quantum-corrected lattice constant of solids, 
i.e. that of computing the harmonic ZPE for a coarse BZ 
sampling ~ preferably only the F-point -, adding it to 
the electronic ground state energy corresponding to the 
equilibrium configuration, and obtaining the lattice con- 
stant as the minimum of the sum. As we shall see now. 
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FIG. 2: (Color online) Contribution of the harmonic ZPE as 
a function the lattice parameter for a 4-atom supercell and 
7=304.97 a.u. Since the ZPE increases as a decreases (red, 
dot-dashed line), the QHA lattice parameter (vertical dashed 
line) corresponding to the minimum of the QHA energy curve 
(blue, dashed), is expanded. 
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FIG. 3: (Color online) Difference in equilibrium lattice con- 
stant with respect to the QHA value, as a function 7 in a 
4-atom ceck 11. The black solid curve is for the ANHA, the 
red (dot-dashed) for the VSCF, and the blue (long-dashed) 
for the exact calculation. 



however, this procedure is not necessarily useful for very 
anharmonic systems such as hydrogen-bonded crystals. 



B. Model hydrogen-bonded chain 

As an extension to highly anharmonic systems, we con- 
sidered a model for a hydrogen-bonded diatomic chain 
where the light atoms are generally hydrogens, and the 
heavy atoms can be oxygen or fluorine. For the sake of 
simplicity, and to be consistent with the following Sec- 



tion, we call them H and F, respectively. For the F-H 
potential we retain a Morse- type of interaction (Eq. [9]) . 
For sufflciently large F-F distances the H atoms feel at- 
tracted to one or the other F-atom, rather than being 
shared between them. Therefore, the H-atoms in the 
chain are subject to a double-well potential. This poten- 
tial has an important property for what concerns isotope 
effects: the barrier between the two wells becomes lower 
when the two neighbouring F-atoms approach each other. 
This is an essential ingredient, which allows for an eas- 
ier migration of the protons upon compression. The F 
displacements (pi) are measured from the given lattice 
parameter. On the other hand, the H displacements (ui) 
are measured from the midpoint between the equilibrium 
positions of the F-atoms. Displacements pi and Ui are 
continuous variables associated to discrete lattice sites i. 
With this, the double Morse contribution to the PES is 

N 2 

£fh = J2 ^ (l-e-''(t+«--f'-'^=)) 
1=1 

N-l 2 

+ Y, D (^i_e-Kt-"'+P'+i-'-<=)) 

i=l 

+ D (^l-e-^(t+""-Pi-'-'=))^ , (10) 

where A'' is the number of sites and a is the lattice pa- 
rameter. The potential parameters are the well depth D 
and the position of the minima re, measured from the 
F-atoms. The effects of periodic boundary conditions 
(PBC) are included in the last term on the RHS. 

When the H-atom sits off-center it forms a neutral unit 
with the F-atom, and another H-atom on the opposite 
side of the F is not welcome due to electronic closed-shell 
repulsion effects. A similar effect occurs in ice, where 
oxygens make two strong covalent bonds with H-atoms 
to form the water molecules, and two weaker hydrogen 
bonds with neighboring molecules. This, however, is not 
taken into account by the Morse potential. To intro- 
duce this feature, usually known as Pauling's "ice-rules" , 
it is necessary to include an interaction between the H- 
atoms that discourages them from approaching the same 
F-atom. There are several possible ways of doing this. A 
popular choice is to include a harmonic interaction be- 
tween H-atoms. The H-H contribution to the PES is 

N-l 

£hH = ^ kffH {u,+i - Uif + kHH {ui - Un)'^ , (11) 
i=l 

where the last term takes care of PBC. This potential 
encourages adjacent H-atoms to move the same amount 
in the same direction, thus promoting the ice rules. The 
larger the constant knH, the more effectively these are 
respected. A similar term is included for the F-atoms 

N~l 

£ff = ^ kpF {Pi+1 - Pif + kFF {pi - PNf , (12) 
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TABLE I: Parameters used for the mode hydrogen-bonded 
hnear chain. AU quantities expressed in a.u. 



D 


b 


kpp 




D, 


bi 


ai 


0.0171 


4.1276 


1.7763 0.02351 


0.06413 


0.80848 


0.3 


4.4 



with the corresponding PBC. This potential favors the 
F-atoms to move in the same direction so that two F- 
atoms cannot approach the same H at once. This model 
has been taken from a work by Yanovitskii et al,^^ who 
studied the phase diagram in the self-consistent harmonic 
approximation for a 0-H- • -O linear chain. The only dif- 
ference with respect to Yanovitskii's model is the value of 
the harmonic constant kHH (0.06413 instead of 0.01413 
a.u.). This modification was necessary because the orig- 
inal value did not enforce the ice rules strongly enough, 
thus leading to a potential with quite unrealistic features 
such as unstable acoustic phonon branches, as shown in 
a preliminary version of this worki^ 

This model assumes that the H-H and the F-F interac- 
tions are independent of the lattice parameter a. A mild 
dependence on the lattice constant enters through the 
H-F interaction, but it is quite unrealistic. In order to 
render quantum effects more realistic we have introduced 
an additional Morse-type potential 



N 

E- 

1=1 



D, 1 



(13) 



The parameters used in the present work were inspired in 
the calculations presented in the following Section, and 
are summarized in Table IH With this choice of param- 
eters, the classical lattice constant is ad = 4.368 Bohr. 
In Fig. |3]we show a schematic picture of the interactions 
involved in this model. 




FIG. 4: (Color online) Schematic view of the interactions in 
a model for a hydrogen-bonded chain. 

Figure [5] shows the equilibrium lattice parameter as a 
function of the inverse number of cells. The behavior 
of the QHA with the number of cells is similar to the 
monoatomic chain (black circles, solid line). Notice that 
here each cell contains 2 atoms, so that a well-converged 
value of QQHA requires at least four cells, i.e. 8 atoms. 
However, three cells are already quite well converged, and 
two cells produce a very reasonable value. Therefore, 



since this is an easily affordable size, in this paper we 
present results mostly obtained with two cells. 

Notice that now aquA decreases with increasing N . In 
fact, at variance with the covalently bonded chain, the 
QHA converged value is 0.025 Bohr smaller than the clas- 
sical value. This is a characteristic feature of H-bonded 
systems, where the stretching frequency that contributes 
to the ZPE initially decreases upon compression because 
the H-bonds weaken. This occurs until the bond becomes 
symmetric, and only then the frequencies (and the ZPE) 
start to increase, as will be shown in Section fill CI 

Next, the normal modes determined at the energy min- 
imum were used to compute the ZPE in the anharmonic 
approximation (ANHA, blue squares, dot-dashed line). 
The converged value of the lattice constant is essentially 
indistinguishable from aquA ■ Its behavior with N , how- 
ever, is non-monotonic. This can be explained through 
the fact that, for this linear chain, there is only one vibra- 
tional coordinate that exhibits a double-well potential. 
As the number of cells is increased, the contribution of 
this double-well to the total energy becomes increasingly 
unimportant. 



4.38 
4.37 
4.36 



b 4.34 

m 

^ 4.33 



4.32 
4.31 
4.30 
4.29 








0.4 0.6 



FIG. 5: (Color online) Quantum-corrected lattice constant for 
a model H-bonded chain as a function of the inverse number of 
cells. Contrary to covalently-bonded systems, the quantum- 
mechanical effect is to decrease a from its classical value (hor- 
izontal long-dashed line on top). The green (light grey, solid) 
line with diamonds represents the VSCF results, while the red 
(short-dashed) line with triangles are exact results. The blue 
(dot-dashed) hne with squares is the ANHA, and the black 
(solid) line with circles is the QHA. 

In order to evaluate the effect of mode-coupling anhar- 
monicities, we carried out VSCF calculations by map- 
ping the PES along the normal modes corresponding to 
the minimum energy configuration, i.e. with all the H- 
atoms off-center (green diamonds, solid line at the bot- 
tom) . The exact calculation used the PES mapped along 
the normal modes corresponding to the saddle-point con- 
figuration, i.e. with all the H-atoms in the center of the 
bonds (red triangles, dashed line). At variance with the 
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VSCF, the choice of vibrational coordinates is irrelevant 
because, unlike the case of molecular systems, the vibra- 
tional subspace in crystals does not depend on the refer- 
ence geometry.™ In the VSCF approximation we went up 
to a 6-atom supercell (A^=3), i.e. a 5-dimensional VSCF 
problem. The good agreement of the VSCF with respect 
to exact results for N = 2, together with the trend with 
N exhibited by the VSCF, suggest that a large fraction 
of the quantum effect is already captured in the two- 
cell calculation. Nevertheless, this statement must be 
taken with caution. Calculations on larger supercells are 
needed to verify this, especially in view of the unusual 
behavior of the anharmonic approximation ANHA. 

Another conclusion is that the QHA underestimates 
the lattice contraction due to quantum effects, remain- 
ing at about half the required value. The next level of 
difficulty and computational cost is the ANHA, but its 
non-monotonic behavior with system size makes it rather 
unsafe as a method to include anharmonicity. In partic- 
ular, the seemingly good results obtained for two cells 
appear to be fortuitous. The VSCF, however, appears to 
be quite a good and robust method that reproduces exact 
results to an excellent accuracy. This indicates that the 
correlation between modes is a relatively minor effect, so 
that a mean-field approximation where each mode feels 
the other modes on average, is sufficient for this class of 
problems. 



C. First- principles hydrogen-bonded chain 

As a realistic application we studied an F-H hydrogen- 
bonded chain by means of first-principles calculations. 
Actual F-H chains develop a zig-zag structure.— How- 
ever, since our purpose here is not to solve the problem 
of F-H crystals, but to understand isotope effects when 
a realistic PES is used, we disregard this geometrical as- 
pect and consider straight chains and the motion of the 
atoms only along the axis of the chain. First-principles 
calculations have been carried out using a combination 
of codes. For phonon calculations we used the pseudopo- 
tential plane wave code Quantum-ESPRESSO.^^ The en- 
ergy curves and lattice constants were calculated with 
SIESTA^ which also uses pseudopotentials but the wave 
functions are expanded in a localized basis set of pseudo- 
atomic orbitals. All calculations were carried out within 
the PBE generalized gradient approximation^^ to density 
functional theory (DFT). 

In Fig. [S] (upper panel) we present the phonon dis- 
persion relations along the direction of the chain, calcu- 
lated in the symmetric configuration with the H-atoms 
in their centered position, from a compressed lattice con- 
stant of 4 Bohr up to an expanded value of 5.2 Bohr. For 
small a there is no double-well, the H-atoms are stable in 
their centered positions and the optical phonon branch 
for this reference configuration is stable. At a w 4.22 
Bohr the double- well emerges and the H-atoms prefer an 
off-centered configuration. If we insist on calculating the 
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FIG. 6: (Color online) First-principles phonon dispersions 
for a linear F-H chain at various lattice parameters. Upper 
panel: H-atoms in centered positions. Red (dashed) lines with 
open squares represent acoustic branches, while black solid 
lines with filled circles correspond to optical branches. Lower 
panel: H-atoms in their stable position. Blue (dashed) lines 
with open squares correspond to stable centered H-atoms, 
while black solid lines with filled circles indicate stable off- 
centered H-atoms. Acoustic branches are omitted for clarity. 
Values of lattice constants in Bohr are indicated in the figure. 



phonon dispersions for the centered configuration, the 
optical branch develops an instability at zone-center. In 
fact, there is a whole portion of the optical branch that 
is unstable. The extent of the instability region increases 
with increasing lattice constant, to the point of making 
the optical mode unstable throughout the whole Brillouin 
zone when the chain is stretched to 5.2 Bohr. 

Further insight is obtained by calculating the phonon 
dispersion relations for the stable configuration. These 
are presented in the lower panel of Fig. [SI where it can 
be observed that the zone-center optical phonon becomes 
soft at a « 4.22 Bohr, and after the H-atoms have cen- 
tered (blue lines and open squares) its frequency starts 
rising again quite steeply. An interesting observation is 
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that the curvature of the optical phonon branch is the 
opposite of the usual picture where the frequency de- 
creases from zone center towards zone boundary. This 
is a consequence of the type of bonding, and common 
to hydrogen-bonded systems. In the zone-center optical 
phonon all unit cells move in phase, with one H-atom 
approaching the F-atom while the other moves away. At 
zone boundary, consecutive unit cells move out of phase 
so that half of the F-atoms are approached by two H- 
atoms rather than one. This is in clear violation of the 
ice-rules, which penalize this type of motion with an in- 
creased energy, thus explaining why the optical mode is 
harder at zone-boundary. 

The harmonic frequencies can now be used to analyze 
the behaviour of the QHA as a function of the lattice 
constant. In Fig. [7] (upper panel) we show the ZPE (red, 
dot-dashed line), the classical PES (black, solid), and 
the quantum-corrected PES (blue, dashed) as a func- 
tion of the lattice parameter for iV = 2 (4 atoms). It 
is interesting to analyze the behavior of the QHA when 
the system is highly compressed. It is apparent that the 
ZPE exhibits a non-analytic behavior (a cusp) at a « 4.2 
Bohr, where the H- atoms center and the type of bonding 
changes from hydrogen-bonding to covalent. As a con- 
sequence, the QHA curve exhibits two minima, one of 
which is spurious. Since this characteristic is associated 
to the stretching mode, and the remaining modes are 
not severely affected by compression, similar behaviour 
is observed for a larger number of cells. This feature rep- 
resents a physically incorrect picture, and thus a clear 
limitation of the QHA in the description of hydrogen- 
bonded systems at high pressures In fact, this non- 
analytic behaviour completely disappears when the prob- 
lem is solved exactly (lower panel). Notice the substan- 
tial difference (0.07 Bohr) between the QHA and exact 
lattice constants (dashed vertical lines). 

Similarly to the model presented in Section fill Bl we 
have examined the dependence of quantum nuclear ef- 
fects on the number of cells, and for the various approx- 
imations. Results are summarized in table HTl below: 

TABLE II: Equilibrium lattice parameters (in Bohr) at vari- 
ous levels of approximation. 

No. cells Classical QHA ANHA VSCF Exact 



4.560 
4.560 



4.195 4.381 
4.462 4.387 



4.381 4.381 
4.390 4.392 



For a single cell there is only one vibrational coordi- 
nate associated to the F-H stretching, which is subject 
to a double-well potential. Since there is only one mode, 
both ANHA and VSCF approximations are equivalent to 
the exact solution. The lattice constant obtained in the 
QHA, however, is anomalously small. By inspection it 
is clear that this value corresponds to the spurious min- 
imum mentioned above, which for one cell appears to 
be lower in energy than the correct minimum. Figure 
[5] shows the quantum-corrected energy as a function of 
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FIG. 7: (Color online) Upper panel: Harmonic ZPE (red, dot- 
dashed line) and QHA corrected energy curve (blue, dashed 
line) . Lower panel: exact ZPE and quantum-corrected energy 
as a function of lattice parameter. Color and line coding as 
in upper panel. The classical energy curve (black, solid line) 
is shown in both panels for comparison. Note the large dif- 
ference between QHA and exact equilibrium lattice constant 
(vertical dashed lines). Arrows indicate the classical lattice 
constant. 



lattice parameter for two cells (i.e. two F-H units, three 
vibrational modes). 

At variance with a single cell, the QHA (red, dot- 
dashed line) provides a better estimate of the effect. 
Using the vibrational coordinates obtained in the sta- 
ble configuration, and already used for the QHA, we in- 
troduced intra-mode anharmonicities through the ANHA 
(black, long-dashed), and mode-coupling anharmonicities 
via the VSCF approach (green, dashed - light grey in 
print). The latter is in very good agreement with the 
exact result (blue, solid), especially for small and large 
values of the lattice parameter. 

In order to assess the quality of the various approxi- 
mations, we mapped the PES along the zone-center op- 
tical mode for different values of the lattice parameter 
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FIG. 8: (Color online) Ground state energy as a function 
of lattice parameter for two unit cells at various levels of ap- 
proximation: The blue (solid) line is for the exact calculation, 
the green (dashed, light grey) is the VSCF, the black (long- 
dashed) is the ANHA, and the red (dot-dashed) the QHA. 



(Fig. [5]). For small values of a the H-atoms are centered 
and subject to a single- well, anharmonic potential (black, 
solid thin line in Fig. [9]). Therefore, at high compression 
the QHA is a good approximation, but closer to the de- 
centering point the potential becomes more anharmonic, 
thus compromising its quality. This situation is notably 
improved by introducing intra-mode anharmonicities in 
the approximation of non- interacting modes (ANHA). 
The VSCF approximation reproduces the exact results 
to an excellent extent, thus indicating that the effect of 
correlation between modes is very small. A similar trend 
is observed for large values of the lattice parameter. In 
this case, the barrier in the double-well potential is so 
high that the overlap between the two degenerate states 
(left and right) is very small. Again, the ANHA pro- 
vides an improved estimate of the energy with respect to 
the QHA, and the VSCF reproduces quite well the exact 
results, thus indicating that correlation between modes 
remains negligible. We conclude that in the two limiting 
cases the ANHA is a rather decent approach, and the 
VSCF reproduces very closely the exact results. 

The situation is somewhat different at intermediate 
values of the lattice parameter, where the barrier is low 
and the effect of anharmonicity is larger. Here the QHA 
and ANHA are not very satisfactory, and have to be im- 
proved by introducing mode-coupling anharmonicities at 
the mean-field level, as in the VSCF scheme. It is also in 
this region where correlations become more important. 

Although the energy curves are shifted, the position 
of the minimum estimated by the ANHA is very good. 
According to the results presented in Section IIII Bl for 
the model hydrogen-bonded chain, we believe this may 
be a fortuitous coincidence for this particular system size. 
The QHA, however, is evidently not sufficiently accurate. 
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FIG. 9: (Color online) First-principles PES along the zone- 
center optical mode for different values of the lattice param- 
eter. The thin black curve represents the most compressed 
situation at a = 4.157 Bohr. The red curve (dashed) corre- 
sponds to a = 4.346 Bohr, the green (long-dashed, light grey) 
to a = 4.535 Bohr (close to the equilibrium value), the blue 
(dot-dashed) to a = 4.724 Bohr, and the magenta (thick solid, 
dark grey) curve to an expanded lattice constant of a = 4.913 
Bohr. Energies are reported for a supercell containing two 
unit cells. 



D. Choice of vibrational coordinates in the VSCF 

The quality of the VSCF approximation depends on 
the choice of vibrational coordinates. The more uncor- 
related they are, the better the approximation. Normal 
modes are a good starting point because in the harmonic 
limit they are completely uncorrelated. However, as the 
amplitude of the displacements along some soft modes 
increases (due to quantum delocalization) , the normal 
modes begin to couple. We now analyze the question of 
which vibrational coordinates are optimal for the VSCF 
approximation, in the case of the 4-atom cell. The three 
modes are the zone-center optical phonon, which exhibits 
a double well, and the zone-boundary acoustic and opti- 
cal phonons, which are practically harmonic. In princi- 
ple, any linear combination of these three modes is possi- 
ble. Nevertheless, since the modes at different k-vectors 
do not mix, we fix the double- well mode and optimize the 
choice of modes in the two-dimensional zone-boundary 
subspacei^ 

We studied the quality of the VSCF for two different 
choices: the normal modes calculated at the minimum of 
the PES, and those calculated at the saddle point config- 
uration, with the H-atoms centered. The VSCF energies 
are reported in Fig. [10] together with the exact energy. 
According to Fig. [51 the two minima in the double-well 
get closer and eventually merge, as the lattice parameter 
is reduced. Therefore, for values of a where the double- 
well has disappeared, there is no saddle point, and the 
normal modes at the minimum are an excellent choice 
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for the VSCF. Similarly, at large a the barrier separating 
the two minima is high, and again the wave function is 
quite localized around the two minima. 




a (Bohr) 

FIG. 10: (Color online) Comparison of the VSCF calculations 
for two different choices of vibrational coordinates. The red 
(dashed) line is for the saddle-point modes, while the black 
(dot-dashed) line is for the modes calculated at the classical 
equilibrium configuration. The blue, solid line represents the 
exact energy. 



Therefore, the vibrational coordinates calculated at 
the minimum of the PES are generally better than the 
saddle-point modes for solving the VSCF problem. At in- 
termediate values of a none of the two choices is clearly 
superior to the other. We have tried to optimize the 
choice of modes within this subspace by minimizing the 
ground state VSCF energy, but the energy gain was al- 
ways of the same order of magnitude of the difference 
between minimum and saddle-point modes. Therefore, 
this appears to be the limit of the VSCF. Any further 
improvement requires the introduction of correlation be- 
tween modes. The present situation is reminiscent of 
static correlation cases often encountered in electronic 
structure calculations, where a single Slater determinant 
(an uncorrelated electronic configuration) is insufficient 
to represent the ground state. These require a multi- 
determinantal wave function as in multi-reference meth- 
ods like CASSCF. Here, a linear combination of two 
Hartree products, each one corresponding to one of the 
two equivalent minima, constitutes an improved wave 
function that correctly respects the symmetry of the 
problem. It is well-known that static correlation cannot 
be recovered by perturbative methods based on a single 
referencei^ This is also a problem for truncated CI ex- 
pansions, while fuU-CI would require a large basis set. 
This is why, in the field of molecular spectroscopy, vari- 
ous groups have developed multi-configuration methods 
such as MCTDH^ 



E. Isotope effects 

Figure [11] shows the effect of modifying the nuclear 
masses on the lattice parameter for the case of two F- 
H cells. In order to evaluate isotope effects, both the 
H and F masses were multiplied by a scale factor s, 
so that the mode eigenvectors, and thus the vibrational 
coordinates, were not modified. We have then solved 
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FIG. 11: (Color online) Energy vs. lattice constant for vari- 
ous mass-scaling factors. The black (solid) line is the classical 
result, blue (dot-dashed) for a — 10, green (long-dashed, light 
grey) for s = 5, violet (dotted) for s = 2, and red (short- 
dashed, dark grey) for s = 1. Vertical lines indicate the equi- 
librium lattice constants, which are then reported in the inset 
as a function of 1/s. The HF chain PES has been calculated 
at the DFT-PBE level, using the SIESTA code. Actual HF 
chains assume a zig-zag form, but here we considered linear 
chains. 

the three-dimensional Schrodinger equation in the vibra- 
tional space for various lattice constants, and obtained 
the quantum corrected energy curves for the true H and 
F masses (red, short-dashed), when these masses are mul- 
tiplied by two (violet, dotted), five (green, long-dashed), 
and ten (blue, dot-dashed). For comparison we also show 
the classical energy curve, corresponding to the limit of 
infinite masses. As expected, the energy increases with 
decreasing mass, while the lattice constant decreases, as 
illustrated in the inset. This is the essence of the geomet- 
ric or Ubbelohde effect in H-bonded systemsj^ Lighter 
particles are more delocalized and exhibit a larger prob- 
ability in the region of the barrier. According to Fig. 
[HI this favors more compressed bonds, thus translating 
into an effective attraction between the two neighbor- 
ing F-atoms (or oxygens in the more common 0-H- • -O 
hydrogen bonds). This attraction is reflected in an en- 
hanced cohesion, and thus a smaller lattice constantj^^ 
which is a common feature of a large family of H-bonded 
crystals. To illustrate this concept of enhanced delocal- 
ization, we show in Fig. [T^j a two-dimensional cut of the 
wave function for the true masses and those multiplied 



11 



by five, for a lattice constant a=4.44 Bohr. 

0.12 r 




FIG. 12: Two-dimensional cut of the ground state wave func- 
tion in the subspace of the two optical modes (2 and (3. The 
upper panel is for the true masses and the lower panel for 
particles five times heavier. Notice the enhanced localization 
of the latter. 

We can now use these wave functions to calculate the 
quantum expectation value of the F-H distance as a func- 
tion of the mass of the particles. This is a very important 
quantity, because it is what is obtained in neutron and 
X-ray diffraction experiments. This distance is usually 
quite different from that obtained classically for the same 
lattice parameter, which corresponds to the minimum of 
the double- well PES. In fact, this is probably the source 
of many inconsistencies between experimental structures 
and those determined via first-principles calculations.^^ 
In Table IIIII we report the exact F-H distance as a func- 
tion of the mass scaling factor for a lattice constant of 
4.44 Bohr. As expected, the F-H distance decreases when 

TABLE III: Equilibrium F-H distance (in Bohr) as a function 
of the mass scaling factor s for a lattice constnat of 4.44 Bohr. 

s 1 2 5 10 (X) 
doH 2.08 2.06 2.02 1.99 1.85 



the masses become heavier. The classical F-H distance 
increases by a substantial 10% when the true H and F 
masses are used. In contrast, at fixed lattice constant the 
QHA values are insensitive to s. This is because quantum 



and classical harmonic oscillators have the same equilib- 
rium position, and thus the quantum harmonic approxi- 
mation does not modify the internal geometry. More gen- 
erally, the F-H distance in the QHA corresponds to the 
classical value obtained at the s-dcpendent equilibrium 
lattice constant, which is reported in the inset to Fig. [Tl] 
(red symbols). Therefore, the influence of mass-scaling is 
limited to the volume effect, which is significantly smaller 
than the full quantum effect. Firstly, the dependence of 
the QHA lattice constant with s is milder than the ex- 
act one. Secondly, as shown in Fig. [51 variations of this 
magnitude lead to small changes in the location of the 
minima of the double-well. As a consequence, the F-H 
distance is severely underestimated in the QHA, thus pre- 
cluding its use in the description of the internal geometry 
of H-bonded systems. 

In the present case, the isotope effect upon doubling 
of all the masses, which is closely related to deutera- 
tion, does not entail a significant change in the F-H dis- 
tance. In a first instance one could think that this is 
because it was calculated at fixed lattice constant (see 
table HIT) ) . In effect, this would be consistent with theo- 
retical calculations^- and experiments under pressurci^ 
When sufficient pressure is applied to deuterated com- 
pounds to reproduce the lattice parameters of its proto- 
nated analogue, the internal geometries turn out to be 
quite similar. Within this context, one of us has shown 
that most of the isotope effect arises from a self-consistent 
interplay between wave function localization, internal ge- 
ometry and lattice parameters. Therefore, to assess the 
full extent of the isotope effect one would have to compare 
the internal geometries at the corresponding equilibrium 
lattice constants. In the present case we obtained F-H 
distances of 2.12 and 2.10 Bohr, respectively. Interest- 
ingly, the isotope effect remains as small as before. This 
is not inconsistent with experimental data, though. In 
fact, the magnitude of the isotope effect on inter-atomic 
distances in hydrogen-bonded systems ranges from al- 
most insignificant values for some systems, to large differ- 
ences as m KDP, where the 0-H distance of 2.45 A rises 
to 2.52 A in its deuterated analogue DKDP.^® This de- 
pends on several factors, mainly the size of the coupling 
between internal coordinates and strain and the shape 
of the double-well (barrier height and distance between 
minima) in the region around the equilibrium lattice con- 
stant. 



IV. CONCLUSIONS AND OUTLOOK 

We have analyzed the influence of quantum nuclear ef- 
fects on the structural properties of solids. To this end 
we have implemented a simple methodology of mapping 
the multidimensional PES in the space of normal modes, 
and then solving the resulting vibrational Schrodinger 
equation. This has been done exactly when possible, and 
also using a number of computationally tractable approx- 
imations that range from the vibrational self-consistent 
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field (VSCF) method and its variants, down to the quasi- 
harmonic approximation (QHA). 

We have studied the behaviour of energies and lat- 
tice constants in an anharmonic monoatomic chain as a 
function of the size of the supercell, and showed that re- 
sults for four cells are probably already quite close to 
convergence. Here, the QHA is an excellent approxima- 
tion unless the masses become exceedingly small, thus 
justifying the customary approach followed in solid-state 
physics. We have then analyzed the behaviour of a model 
hydrogen-bonded chain, and observed that the coupling 
between modes can be important. Finally, we consid- 
ered a realistic H-bonded chain where the PES was de- 
termined from first-principles calculations. Here we ob- 
tained the quantum-corrected lattice constant and inter- 
nal structural parameters (distances), and showed how 
isotope effects arise in these systems. We have also 
analyzed the various approximations and showed that 
the QHA is insufficient, while the independent-mode an- 
harmonic approximation ANHA appears to introduce 
an important improvement at the energetic and struc- 
tural level. Nevertheless, results obtained for the model 
hydrogen-bonded chain suggest that this may be acciden- 
tal and due to size effects. VSCF results are very close 
to exact ones, thus suggesting that correlation between 
modes is a minor effect. Optimizing the vibrational coor- 
dinates in the VSCF method does not constitute a signif- 
icant advantage; further improvement should arise from 
correlated methods such as vibrational CI or MCTDH. 
In any case, the VSCF approach using the modes calcu- 
lated at the classical equilibrium geometry appears to be 
a very good approximation. 

Although technically the extension to higher dimen- 
sionalities is straightforward, there are some features that 
are not present in one-dimensional systems. For exam- 
ple, the tunneling mode can be strongly coupled to other 
modes^ so that, for large amplitude of motion, correla- 
tion between them may become important. In this case 
we enter the realm of multi-dimensional tunneling where, 
if we insist on describing the proton motion in terms of 
a single, effective tunneling coordinate, this latter tends 
to be curvilinear. Therefore, it is quite likely to have to 
adopt schemes that combine the exact solution for one 
or a few subspaces of strongly coupled modes, with a 
VSCF representation or a lower-level scheme such as the 
QHA for the remaining modesi^ It remains the problem 
of identifying which are those subspaces, though. 

There are many aspects of this approach that can be 
improved, especially for what concerns simplified meth- 
ods that can be used for larger systems and higher dimen- 
sionality. Within the VSCF scheme, the limiting factor 



is the mapping of the PES and the multi-dimensional 
integration. The latter can be efficiently dealt with by 
factorizing the PES into a sum of products of modesi^ 
This allows for the computation of VSCF potentials as 
products of inexpensive one-dimensional integrals. The 
problem of fitting the PES remains open, although there 
are general strategies for PES fitting based on interpola- 
tion methods, ■^^ or using principles of multivariate analy- 
sis such as the high-dimensional model representation of 
Rabitz and Ali§j^ An improvement to the ANHA that 
approximately recovers the interaction between modes at 
the mean-field level consists of using the quantum expec- 
tation value of the geometry as a reference, rather than 
the classical one. This can also be seen as a simplification 
of the VSCF approximation where the vibrational wave 
functions in the integrals are replaced by completely lo- 
calized delta-functions. This is a good approximation 
unless one or more approximated wave functions exhibit 
two peaks corresponding to a double-well potential. 

Another subtle issue in hydrogen-bonded systems is 
that, although DFT-GGA approaches reproduce struc- 
tural properties quite well, proton transfer barriers 
appear to be exceedingly small, in some cases even 
disappearing altogetherj^^, Therefore, other electronic 
structure methods such as hybrid Hartree-Fock-DFT— 
or correlated quantum chemical methods have to be 
explored 1^ 

A straightforward extension of this methodology is the 
calculation of vibrational excitations. Once the PES is 
known, excited states can be easily calculated and used 
to compute thermodynamic quantities beyond the quasi- 
harmonic approximation. General VSCF algorithms use- 
ful to tackle mild anharmonicities in molecular systems 
have been implemented many years ago,^^ and more re- 
cently in conjunction with an ab initio description of the 
PES."^^ In fact, a VSCF option is available in some elec- 
tronic structure codes like GAMESS.^ The mam use 
of this capability is at the spectroscopic and thermo- 
chemical, rather than structural level. An extension of 
this methodology to crystalline systems is straightfor- 
ward and perfectly viable, although to the best of our 
knowledge this has not been implemented so far. Never- 
theless, the treatment of highly anharmonic systems such 
as those involving double- wells remains a challenge. We 
expect the present paper to be a relevant contribution in 
this direction. 
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